function sim = solve_bf(sim, eq_SS, param, glob, options)

% Solves "backwards" and iterates "forwards" on one iteration of the MIT
% shock

% start at the end and solve backwards
v2    = reshape(eq_SS.v.v2, glob.n(1), glob.n(2)); % assume that the expected value is in SS at date T

tic() 

for t = options.T:-1:1
   param.M = sim.M(t);
    
    param.C = sim.C(t);
    param.D = sim.D(t);
    param.w = sim.W(t); % HH FOC
    param.A = sim.A(t); 
    
    
    c = reshape(v2, glob.n(1), glob.n(2));
    c = [c; c];

    
    % First solve the incumbent's value function
    v     = solve_valfunc(c,glob.s,param,glob,options);
    
    v2    = reshape(v.v2, glob.n(1), glob.n(2));

    % Solve on a finer grid to get policy functions
    tmp      = glob.PnK;
    glob.PnK = glob.PnKf;
    vf       = solve_valfunc(c,glob.sf,param,glob,options);
    glob.PnK = tmp;
 
    if t==1
        disp('pausing')
    end
    
    % Record optimal policy
    sim.y(:,t)  = vf.y; 
    sim.l(:,t)  = vf.l;
    sim.p(:,t)  = vf.p;  
    
    sim.v(:,t)  = vf.v1;
        
    if (t == 1)
        disp('hi')
    end
    
end
toc()

%% Now iterate forwards on policies
tic()

fspaceerg   = fundef({'spli',glob.bgridf,0,1});


for t = 1:options.T
    
    param.M = sim.M(t);

    % Distribution from last period
   if t > 1
       L               = sim.mu(:,t-1);
       
       lp = sim.l(:,t-1);       
   else
        L               = eq_SS.L;
        
        lp = eq_SS.l;
   end
      
   
    % set up transition matrix
    lp          = min(lp, glob.maxb);
    lp          = max(lp, glob.minb);
      
    Qb              = funbas(fspaceerg, lp);
    Q               = dprod(glob.Qef,Qb);
        
    sim.mu(:,t) = Q'*L; 

    eq_mit.L          = sim.mu(:,t);
    eq_mit.yf         = sim.y(:,t);
    
    aggs              = find_agg_mit(eq_mit, param, glob, options);
    
    sim.C(t)          = aggs.C;
    sim.D(t)          = aggs.D;

    L                 = sum(eq_mit.L.*sim.l(:,t)); % this is actually the previous period's labor supply - FIXED WLG 7/30                    
    sim.L(t)          = L;
    sim.W(t)          = L^param.nu*param.psi;
    
end
toc()

% Labor       =  sim.l;
% sim.L       =  sum(sim.mu.*Labor);
% sim.W       =  sim.L.^param.nu*param.psi;


end